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I. INTRODUCTION 



Colloidal suspensions are of fundamental interest for various applications. One of the ba- 
sic problems of colloidal science is how to stabilize a lyophobic colloidal suspension against 
flocculation and precipitation. A common approach is to synthesize particles with acidic 
or basic charged groups on the surface 1 -^. When placed in a polar medium such as water, 
these groups become ionized and the particles acquire a net charge. Repulsion between 
like-charged colloidal particles then prevents them from coming into a close contact where 
the short-range van der Waals forces become important. Addition of electrolyte to colloidal 
suspension leads to screening of the Coulomb repulsion 3 . At critical coagulation concen- 
tration (CCC), the repulsive energy barrier disappears and the van der Waals forces drive 
colloidal coagulation and precipitation 4-7 . It is also well known that addition of even very 
small amount of multivalent ions leads to a rapid precipitation. The correlation induced 
attraction between the colloidal particles produced by the multivalent ions is sufficient to 
precipitate colloidal suspensions even without taking into account the van der Waals forces. 
The like-charge attraction has been extensively explored in colloidal and polyelectrolyte lit- 
erature^—. A related phenomenon known as the charge reversal has also attracted a lot of 
attention over the recent years 1 ^— . In this case electrostatic correlations result in a strong 
colloid-counterion association 3 . The counterion condensation can be so significant as to 
reverse the electrophoretic mobility of colloidal particles^ 1 ^. 

Most of the theoretical work on stability of colloidal suspensions and the charge reversal, 
however, neglects the effects of the dielectric discontinuity at the particle/solvent interface. 
In fact, in many colloidal suspensions the static dielectric constant of colloidal particles 
can be 20 to 40 times lower than the static dielectric constant of the surrounding water. 
This means that an ion in the vicinity of colloidal surface will encounter a strong ion-image 
repulsion. This repulsion can significantly affect the effective charge of the colloid-counterion 
complex and thus modify the colloid-colloid interaction potential. The polarization effects, 
however, have been mostly neglected in almost all of the theoretical studies. The reason 
for this is that it is very hard to include the dielectric discontinuities in anything but the 
simplest planar geometry. Thus, even to perform a Monte Carlo simulation that accounts 
for the dielectric discontinuity requires a significant computational effort. Some years ago 
Linse^i proposed to account for the induced charges by treating the low dielectric colloidal 
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particle as if it was an "inverse" conductor. It is well known that if one places a charge near 
a conducting sphere, a surface charge will be induced on the sphere^. The field produced 
by the surface charge is exactly equivalent to the field produced by two point charges of 
opposite sign, one located at the spheres inversion point and another at its center. In the 
case of a conductor, the charge at the inversion point has the opposite sign to the charge 
placed outside the sphere, so that this charge is attracted to the conductor. Linse suggested 
that the low dielectric sphere in water can be treated as an inverse conductor, meaning 
that the same construction should apply to locate the images, but that their sign will be 
the opposite of the images inside the conducting sphere. This, however, is not quite right. 
Because of the complicated boundary conditions (BCs) imposed by the Maxwell equations 
at the dielectric interface, one can not satisfy the BCs with only two pointlike image charges. 
In fact one needs an infinite number of images^ 1 ^. 

Polarization effects have also been explored in ionic liquids—"— and for polyelectrolyte 
adsorption^ - — in a slab geometry. In this paper we will derive the explicit inter-ionic inter- 
action potential which very accurately accounts for the dielectric discontinuity for spherical 
colloidal particles. We compare our results with the Monte Carlo simulations of Messina^, 
who obtained the ion-ion interaction potential as an infinite series in Legendre polynomials. 
In the final part of the paper we will analyze the effect of 1:1 electrolyte on the distribution 
of trivalent and monovalent counterions near the colloidal surface. 

II. METHOD 

We will use a primitive model of a colloidal suspension in which colloidal particle is 
represented by a sphere of radius a and the dielectric constant e c . Water will be modeled 
as a uniform dielectric of permittivity e w . The system is at room temperature, so that the 
Bjerrum length, defined as A# = q 2 /e w kBT, is 7.14A. Consider an a-valent ion of charge 
Q = aq, where q is the proton charge, at position from the center of the colloidal particle. 
The Maxwell equations require the continuity of the tangential component of the electric field 
and the continuity of the normal component of the displacement field across the colloid-water 
interface. It is possible to show^ 1 ^ that this boundary conditions can be satisfied exactly by 
placing an image charge Q' = 'yQa/ri inside the colloid at the inversion point r'j = r^a 2 jr 2 
and a counterimage line-charge A (it), distributed along the line connecting the center of 
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FIG. 1. Illustrative representation of an ion of charge Q near a colloidal particle and the induced 
image, Q', and counterimage, A(it), charges. 

colloid with the inversion point r-, with line-charge density 

A(M) = ~5T" UJ ' (1) 

where 7 = (e to — e c ) / (e w + e c ) and < u < r • is the distance along the line, see Fig. [TJ Note 
that this construction does not change the net charge of the colloidal particle, that is the 
total counterimage charge is —Q'. 

The electrostatic potential produced by the image charge at an arbitrary position r outside 
colloid is 

ip im {r; n) = ^ 9a a2 , (2) 

i 

and the electrostatic potential produced by the counterimage line-charge is 

a 2 f l HVt) 
^ Ci (r;ri) = / drj- ^— . (3) 

i 

For r = r i; the integral can be performed exactly in terms of the hypergeometric function 
2-F1. We find the counterimage-ion interaction potential to be 



C'W = 2^1 ( J + ^7 1^ + ;7 3 ) ■ W 
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Although Eqs. E] and H] are exact, they are not very useful for Monte-Carlo or molecu- 
lar dynamics simulations — the integral in Eq. [3] must be done numerically for each new 
configuration of ions, making simulations very slow. However, we can consider a simpli- 
fying approximation. We note that the dielectric constant of a colloidal particle is much 
smaller than the dielectric constant of the surrounding medium. Thus, to leading order 



in e c /e w we can take 7 « 1. In this case the counterimage charge is uniformly distributed, 
\{u) = —Q'/r'i, and the integral in Eq. [3]can be performed exactly, yielding the counterimage 
potential at an arbitrary position r, 

t / x aq ( rri-r-Y; \ 

w ci r: rj = log , 5 

V ; e w a 5 \a 2 - r • n + ^a 4 - 2a 2 (r ■ r ; ) + r 2 rf J KJ 

where the over-bar is used to denote the uniform line-charge approximation. The ion- 
counterimage interaction potential also reduces to a simple equation, 

e,„a " \ r 2 



C (/ W = - log(l-^) • (6) 



III. MONTE CARLO SIMULATIONS 

The simulations are performed inside a spherical Wigner-Seitz (WS) cell of radius R with 
a colloidal particle of charge — Zq placed at the center. The cell also contains N = Z/a 
ct-valent counterions each of diameter d. The electrostatic potential produced at position r 
by an ion located at 1^ is 

0( r ; r i) = —r^- — r + 7" ° tt a , + 7^c»(r; r«) , (7) 

i 

where the first term is the electrostatic potential produced by the ion and the second and 
the third terms are the potentials produced by the image and the counterimage charges, 
respectively. The first two terms of Eq. [7| are exact. In the third term we have used the 
condition of charge neutrality to correct the ion-counterimage interaction from Eq. [5] by 
including a prefactor 7 in front of ip C i{ r ', r i)- This, then, is the Green function for the present 
geometry. The interaction potential between two ions i and j is ag^r^r,). The work 
required to bring all the ions from infinity to their respective positions inside the cell is, 



N ^ 9 N N-1 N 

L 

i=l w * i=l i=l j=i+l 



U = H~ + Y: u '" f + £ E a ^ r J • ( 8 ) 

ton I 



TT seif _ ia 2 q 2 a ag^ 8 ?* (n) 

where U 8el ^ is the interaction energy of the ion i with its image and counterimage charges. 

We use the Eq. [8] in a typical Metropolis algorithm 33 , with 10 5 MC steps for equilibration 
and 10 4 steps for production. We obtain the ionic density profiles dividing the WS cell in 
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FIG. 2. The dashed lines represent the density profiles obtained using the present method and 
the symbols represent the profiles obtained by Messina*^. The parameters of the simulations are: 
e w = 80, e c = 2, d = 3.57A, a = 7.5d, R = AOd, Z = 60, and r c = a + d/2 is the contact distance. 

volumetric bins and counting the average number of particles in each bin for all uncorrelated 
configurations. In Fig. [2] the profiles are compared with the ones obtained by Messina^. 
The agreement between the two simulations is excellent, with a huge gain in computational 
time. To take into account the correct boundary conditions, Messina calculated the inter- 
ionic interaction potential as an infinite series in Legendre polynomials. This method is 
extremely slow, since one needs to calculate hundreds of terms of the infinite series for each 
new configuration in order to obtain a good convergence. To compare the time of processing, 
we perform the simulation of Messina for parameters described in Fig. [2] with trivalent ions. 
Even for a very small number of counterions — twenty counterions — used by Messina, 
expansion in Legendre polynomials is ~ 458 x slower than the method presented in the 
current paper. To speed up the simulations Messina tabulated the counterion-counterion 
interaction potential. Nevertheless, his approach remains at least one order of magnitude 
slower than our Green function method, and is significantly more difficult to extend to larger 
system sizes. 

The approximation of the counterimage charge by a uniform line-charge density, should 
work very well for colloids of low dielectric constants, which are of most practical inter- 
est. However, it is interesting to examine up to what value of e c does this approximation 
remains accurate. Using the exact numerical evaluation of the integral in Eq. [3J and the 
hypergeometric representation of the counterimage-ion interaction potential Eq. HJ we have 
performed the simulations for different values of e c using the exact numerically calculated 
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FIG. 3. Density profiles for various e c . The symbols represent the density profiles obtained using 
the exact counterimage line-charge distribution, while the dashed lines are calculated using the 
approximate uniform counterimage line-charge distribution, Eqs. [5]and[6l The parameters of the 
simulations are the same as in Fig. [2j for a = 3. 

interaction potential and compared the counterion density profiles with the ones obtained in 
simulations with the approximate interaction potentials, Eqs. [5] and [6j In Fig. [3] we show the 
results of these simulations. We see that the approximation works very well up to e c ~ 20, 
which are in the range of most practical interest. 

The simulations using the method developed in the present work are so quick that it is easy 
to study systems which contain mixtures of multivalent and monovalent electrolytes. We 




FIG. 4. Density profiles of trivalent counterions for various concentrations of 1:1 electrolyte. The 
parameters of the simulations are the same as in Fig. [2] for R = 25d, a = 1 and Z = 40. The 3:1 
concentration is 20 mM. 
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next consider a WS cell that contain 3:1 electrolyte at concentration p t and 1:1 electrolyte at 
concentration p m . The number of trivalent counterions inside the system is N t = pt^-(R 3 — 
a 3 ), the number of monovalent counterions is N m = p m ^(R 3 — a 3 ) and the number of 
monovalent coions is AL = 3N t + N m . In Fig. HI the density profiles of 3:1 salt cations are 
presented for various concentration of 1:1 salt. As expected, with increase of the monovalent 
salt concentration, more trivalent cations prefer to be solvated in the bulk of suspension, 
where their electrostatic self-energy is screened most effectively by the other ions^ 1 ^. 
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FIG. 5. Maximum density of trivalent salt counterions as function of concentration of 1:1 elec- 
trolyte. The parameters of the simulations are: e w = 80, e c = 0, d = 4A, a = 30A, R = 70A and 
Z = 90. 

The charge-image repulsion results in density profiles of trivalent ions which have a char- 
acteristic maximum near the colloidal surface. In Fig. we examine the effects of 3:1 and 
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FIG. 6. The density at contact of monovalent counterions for varying concentrations of 3:1 elec- 
trolyte. The parameters of the simulations are the same as in Fig. [5j 
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1:1 electrolyte on the maximum density of trivalent counterions near the colloidal surface. 
Again we see that increasing the concentration of 1:1 electrolyte diminishes the counterion 
condensation — resulting in a smaller counterion density in the vicinity of the colloidal 
surface. More surprising, perhaps, is the behavior of the contact density of the monovalent 
counterions, Fig. |HJ We see that at small concentrations of 3:1 electrolyte the contact den- 
sity varies significantly with the concentration of 1:1 electrolyte. This dependence, however, 
rapidly saturates, so that for 50 mM of 3:1 electrolyte, we no longer see any variation of the 
contact density with the concentration of 1:1 salt. The Fig. [6]shows that with increasing 3:1 
concentration the condensed monovalent counterions are rapidly replaced by the trivalent 
ones. 

IV. CONCLUSIONS 

We have presented a very efficient method for simulating colloidal suspension composed 
of lyophobic colloidal particles of low dielectric constant. The method relies on the exact 
calculation of the Green function for the spherical geometry. The results are in excellent 
agreement with the earlier simulations of Messina^ 2 - — who used expansion in Legendre 
polynomials to account for the dielectric discontinuity at the colloidal surface — with a huge 
gain in the computation time. At the moment, we have only implemented the simulation 
inside a WS cell geometry. In the future, an effort should be made to extend the theory to 
take into account periodic boundary conditions through the use of Ewald summation. 

This work was partially supported by the CNPq, Fapergs, INCT-FCx, and by the US- 
AFOSR under the grant FA9550-09- 1-0283. 
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